Kinetic Models in Biology and Chemistry

Robert B. Nachbar, Irene Pak, and Andrew (Spann^1)

Merck Research Laboratories
126 E. Lincoln Ave.
Rahway, NJ 07065

\(TraditionalForm\`\(\(\ \)\(\^1\)\)\) Current address: Massachusetts Institute of Technology.

Abstract

Dynamically changing systems in biology and chemistry are often modeled as a system of ordinary differential equations (ODEs). Biologists and chemists often represent these systems as a network of entities connected by arrows. For example, the Michaelis-Menten enzyme catalyzed reaction is written thus

({{\(Enz + Sub\( ⇌ \+\(k\_off\)\%\(k\_on\)\)ES\)}, {ES (Overscript[→ , (Subscript[k, p])])Prod}})

The model for absorption and distribution of a drug is often drawn pictorially, as shown on the left, and can also be represented mechanistically, as on the right,

({{, ({{ (Overscript[↦ , dose])(Subscript[drug, plasma])}, {\(drug\_gut\( ⇒ \&\(k\_a\)\)drug\_plasma\)}, {\(drug\_plasma\( ⇔ \+\(k\_21\)\%\(k\_12\)\)drug\_tissue\)}, {(Subscript[drug, plasma]) (Overscript[⇥ , (Subscript[k, e])])}})}})

We will describe the package ChemicalKinetics.m that facilitates the use of these mechanistic representations for constructing the system of ODEs for both mass-action and Michaelis-Menten models. We have also developed two other packages: UnitsTools.m and ODEDerivatives.m. The functions in the former help the user verify that the dimensional units of the parameters of the model are self-consistent (and thereby avoid problems with transferability), and make it easy to use units with numerical functions such as NDSolve, NMinimize, and NonlinearRegress. The functions in the latter package define the rules by which Mathematica can automatically take derivatives with respect to parameters of the system of ODEs so that the user can perform sensitivity analyses and use gradient-based methods for parameter fitting.

Modeling hierarchy

Modeling system requirements

Three packages

ChemicalKinetics.m

UnitsTools.m

ODEDerivatives.m

ChemicalKinetics.m

\[KernelIcon]

<< Chemistry`ChemicalKinetics`

({{Palette, Functions}, {({{, ({{ }, {Reactions and equlibria}, {Diffusion between compartments}, {Drug administration and elimination}, {Compartment notation}, {Concentration  & amount notation}})}}), ({{ResolveModel, Model integrity}, {KineticModel, Translate model to system of ODEs}, {DimensionalQ, Variable and units consistency}})}})

Overall, balanced reaction vs. stepwise model

An apparently simple reaction

2  NO + 2  (Subscript[H, 2]) → (Subscript[N, 2]) + 2  (Subscript[H, 2]) O

may be a series of elementary reactions and equlibria

({{\(2  NO\( ⇌ \+\(k\_2\)\%\(k\_1\)\)\(N\_2\) O\_2\)}, {(Subscript[N, 2]) (Subscript[O, 2]) + 
              (Subscript[H, 2]) (Overscript[→ , (Subscript[k, 3])])(Subscript[N, 2]) O + (Subscript[H, 2]) O}, {(Subscript[N, 2]) O + (Subscript[H, 2]) (Overscript[→ , (Subscript[k, 4])])(Subscript[N, 2]) + (Subscript[H, 2]) O}})

Mass action kinetics

\[KernelIcon]

\(\(model = {2  NO\( ⇌ \+\(k\_2\)\%\(k\_1\)\)N2O2, 
          N2O2 + H2\( → \&\(k\_3\)\)N2O + H2O, 
          N2O + H2\( → \&\(k\_4\)\)N2 + H2O};\)\)
{vars, {odes, ics}} = 
        KineticModel[model, t, 
          InitialConditions -> {[NO][
                  0] == 𝒩(Subscript[𝒪, 0]), [H2][
                  0] == ℋ(Subscript[2, 0])}];

\[KernelIcon]

odes // ToColumns[#, 1] & // 
      Grid[#, RowSpacings -> 3, 
          ColumnAlignments -> "=="] &

\[KernelIcon]

ics // ToColumns[#, 2] & // 
      Grid[#, RowSpacings -> 3, 
          ColumnAlignments -> "=="] &

ResolveModel

Sometimes, either inadvertantly or through an error, the model is separable. That is, it consists of two or more sets of reactions that do not share any reactants or products.

\[KernelIcon]

\(\(model2 = {H2 + N2O2\( → \&\(k\_1\)\)H2O + 
              N2O, \(HCO3\^-\) + \(Na\^+\)\( ⇌ \+\(k\_3\)\%\(k\_2\)\)NaHCO3, 
          2\ NO\( ⇌ \+\(k\_5\)\%\(k\_4\)\)N2O2, 
          H2 + N2O\( → \&\(k\_6\)\)H2O + N2, 
          NaOH\( ⇌ \+\(k\_8\)\%\(k\_7\)\)\(OH\^-\) + \(Na\^+\), 
          CO2 + \(OH\^-\)\( ⇌ \+\(k\_10\)\%\(k\_9\)\)\(HCO3\^-\)};\)\)
{vars2, {odes2, ics2}} = KineticModel[model2, t];

\[KernelIcon]

ResolveModel[model2]

Rate constants

Values for rate constants are inserted into the ODEs with ReplaceAll and a list of Rules. The user can use whatever units are most convenient, although units from the standard add-on package Miscellaneous`Units` are preferred.

\[KernelIcon]

params = {(Subscript[k, 1]) -> 120.  (Second^(-1)), 
          (Subscript[k, 2]) -> 300.  (Second^(-1)), 
          (Subscript[k, 3]) -> 180.  ((((Mole/Liter)))^(-1)) (Second^(-1)), 
          (Subscript[k, 4]) -> 
            4.  ((((Mole/Liter)))^(-1)) 
              (Second^(-1)), 𝒩(Subscript[𝒪, 0]) -> ((0.5  Mole)/(22.  Liter)), ℋ(Subscript[2, 0]) -> ((0.5  Mole)/(22.  Liter))};

DimensionalQ

When units are included with the parameters, we can test for internal consistency.

\[KernelIcon]

DimensionalQ[odes /. params,  ics /. params,  t]

We can see that the units for (Subscript[k, 1]) are for a unimolecular reaction instead of a bimolecular reaction.

\[KernelIcon]

params={(Subscript[k, 1])->120.((((Mole/Liter)))^(-1))(Second^(-1)),(Subscript[k, 2]) -> 300.  (Second^(-1)),(Subscript[k, 3]) -> 
              180.  ((((Mole/Liter)))^(-1)) (Second^(-1)),(Subscript[k, 4]) -> 
              4.  ((((Mole/Liter)))^(-1)) (Second^(-1)),𝒩(Subscript[𝒪, 0]) -> ((0.5  Mole)/(22.  Liter)),ℋ(Subscript[2, 0]) -> ((0.5  Mole)/(22.  Liter))};
DimensionalQ[odes /. params,  
      ics /. params,  t]

UnitsTools.m

Another common problem is mixing units for a given dimension (remember the Mars Climate Orbiter crash in 1999?). To address this problem and others, we have developed tools to facillitate working with units.

\[KernelIcon]

?ShedUnits
?RestoreUnits
?StandardizeUnits

\[KernelIcon]

SetOptions[StandardizeUnits,  StandardAmount -> Mole, 
      StandardTime -> Second]

For example, using both Second and Minute for time:

\[KernelIcon]

{(Subscript[k, 1]) -> 
            120.  ((((Mole/Liter)))^(-1)) (Second^(-1)),(Subscript[k, 2]) -> 300.  (Minute^(-1)),(Subscript[k, 3]) -> 
            180.  ((((Mole/Liter)))^(-1)) (Second^(-1)),(Subscript[k, 4]) -> 
            4.  ((((Mole/Liter)))^(-1)) (Second^(-1)),𝒩(Subscript[𝒪, 0]) -> ((0.5  Mole)/(22.  Liter)),ℋ(Subscript[2, 0]) -> ((0.5  Mole)/(22.  Liter))}//StandardizeUnits

Running a model

\[KernelIcon]

soln = 
        NDSolve[{odes, ics} /. params // ShedUnits, vars, {t, 0, 600}] // 
          First;
plotModel[soln];

Michaelis-Menten enzyme kinetics

The model for an enzyme catalyzed reaction is

({{E, +, S, \( ⇌ \+\(k\_\(-1\)\)\%\(k\_1\)\), ES,  (Overscript[→ , (Subscript[k, p])]), E + P}})

where S is the substrate, E is the enzyme, ES is the enzyme-substrate complex, and P is the product. The ODEs are

\[KernelIcon]

\(\({vars, {odes, ics}} = 
        KineticModel[{S + 𝔼\( ⇌ \+\(k\_\(-1\)\)\%\(k\_1\)\)ES, 
            ES\( → \&\(k\_p\)\)𝔼 + P}, 
          t];\)\)
odes // ToColumns[#, 1] & // 
      Grid[#, RowSpacings -> 3, 
          ColumnAlignments -> "=="] &

By making the steady-state  assumption that ([ES])'[t]==0, conserving mass, and some substitutions of rate constants for equilibrium constants, one can derive the Michaelis-Menten rate for [P][t].

\[KernelIcon]

soln=Solve[{
([S])'[t]==(Subscript[k, -1]) [ES][t] - 
                  (Subscript[k, 1]) [S][t] [𝔼][t],
([P])'[t]==(Subscript[k, p]) [ES][t],

0 == -(Subscript[k, -1]) [ES][
                      t] - (Subscript[k, p]) [ES][t] + 
                  (Subscript[k, 1]) [S][t] [𝔼][t]   ,

[(Subscript[𝔼, T])][
                  t] == [𝔼][t] + [ES][
                    t]   ,
(Subscript[K, S]) == (([S][t] [𝔼][t])/([ES]))[t]   ,
(Subscript[V, max]) == 
                (Subscript[k, p])×[(Subscript[𝔼, T])][t]   },
{([S])'[t],([P])'[t]},
{[ES][
              t], [𝔼][
              t], [(Subscript[𝔼, T])][t], (Subscript[k, 1]), (Subscript[k, -1]), 
            (Subscript[k, p])}]

We have implemented a short-hand notation for this kind of reaction

\[KernelIcon]

model = {S (Overscript[→ , MM[enz, (Subscript[K, S]), 
                  (Subscript[V, max])]])P};
{vars, {odes, ics}} = 
        KineticModel[model,  t, 
          InitialConditions -> {[S][
                  0] == (Subscript[𝒮, 0])}];
odes

Michaelis-Menten inhibitor kinetics

The mechanism for an inhibitor of the enzyme that binds competitively to the same site as the substrate and its rate laws are

({{({{E, +, S,  (Overscript[⇌ , (Subscript[K, S])]), ES,  (Overscript[→ , (Subscript[k, p])]), E + P}, {+,  ,  ,  ,  ,  ,  }, {I,  ,  ,  ,  ,  ,  }, {    ⥮(Subscript[K, I]),  ,  ,  ,  ,  ,  }, {EI,  ,  ,  ,  ,  ,  }}), ({{([S])'[t]->\(-\(\(V\_max\ \([S]\)[t]\)\/\(K\_S + \([S]\)[
                            t] + \(K\_S\ \([𝕀]\)[t]\)\/K\_𝕀\)\)\)}, {([P])'[t]->\(\(V\_max\ \([S]\)[t]\)\/\(K\_S + \([S]\)[
                          t] + \(K\_S\ \([𝕀]\)[t]\)\/K\_𝕀\)\)}, {([𝕀])'[t]->0}})}})

An inhibitor that binds to the enzyme-substrate complex (uncompetitive inhibition) has this mechanism and rate laws

({{({{E, +, S,  (Overscript[⇌ , (Subscript[K, S])]), ES,  (Overscript[→ , (Subscript[k, p])]), E + P}, { ,  ,  ,  , +,  ,  }, { ,  ,  ,  , I,  ,  }, { ,  ,  ,  ,     ⥮(Subscript[K, I]),  ,  }, { ,  ,  ,  , IES,  ,  }}), ({{([S])'[t]->\(-\(\(V\_max\ \([S]\)[
                            t]\)\/\(K\_S\ \((1 + \([S]\)[t]\/K\_S + \(\([S]\)[t]\ \([𝕀]\)[t]\)\/\(K\_S\ K\_𝕀\))\)\)\)\)}, {([P])'[t]->\(\(V\_max\ \([S]\)[
                          t]\)\/\(K\_S\ \((1 + \([S]\)[t]\/K\_S + \(\([S]\)[t]\ \([𝕀]\)[t]\)\/\(K\_S\ K\_𝕀\))\)\)\)}, {([𝕀])'[t]->0}})}})

Noncompetitive inhibition has this mechanism and rate laws

({{({{E, +, S,  (Overscript[⇌ , (Subscript[K, S])]), ES,  (Overscript[→ , (Subscript[k, p])]), E + P}, {+,  ,  ,  , +,  ,  }, {I,  ,  ,  , I,  ,  }, {    ⥮(Subscript[K, I]),  ,  ,  ,     ⥮(Subscript[K, I]),  ,  }, {IE, +, S,  (Overscript[⇌ , (Subscript[K, S])]), IES,  ,  }}), ({{([S])'[t]->\(-\(\(V\_max\ \([S]\)[
                            t]\)\/\(K\_S\ \((1 + \([S]\)[t]\/K\_S + \([𝕀]\)[t]\/K\_𝕀 + \(\([S]\)[t]\ \([𝕀]\)[t]\)\/\(K\_S\ K\_𝕀\))\)\)\)\)}, {([P])'[t]->\(\(V\_max\ \([S]\)[
                          t]\)\/\(K\_S\ \((1 + \([S]\)[t]\/K\_S + \([𝕀]\)[t]\/K\_𝕀 + \(\([S]\)[t]\ \([𝕀]\)[t]\)\/\(K\_S\ K\_𝕀\))\)\)\)}, {([𝕀])'[t]->0}})}})

Which in our short-hand are

\[KernelIcon]

model = {S (Overscript[→ , MM[enz, (Subscript[K, S]), 
                  (Subscript[V, max])]])P, 𝕀 (Overscript[⇥ , MMcomp[enz, 
                  (Subscript[K, 𝕀])]])};
{vars, {odes, ics}} = 
        KineticModel[model,  t, 
          InitialConditions -> {[S][
                  0] == (Subscript[𝒮, 0])}];
odes

\[KernelIcon]

model = {S (Overscript[→ , MM[enz, (Subscript[K, S]), 
                  (Subscript[V, max])]])P, 𝕀 (Overscript[⇥ , MMun[enz, 
                  (Subscript[K, 𝕀])]])};
{vars, {odes, ics}} = 
        KineticModel[model,  t, 
          InitialConditions -> {[S][
                  0] == (Subscript[𝒮, 0])}];
odes

\[KernelIcon]

model = {S (Overscript[→ , MM[enz, (Subscript[K, S]), 
                  (Subscript[V, max])]])P, 𝕀 (Overscript[⇥ , MMnon[enz, 
                  (Subscript[K, 𝕀])]])};
{vars, {odes, ics}} = 
        KineticModel[model,  t, 
          InitialConditions -> {[S][
                  0] == (Subscript[𝒮, 0])}];
odes

General form for Michaelis-Menten inhibitor kinetics

Given a general system of (Subscript[n, s]) substrates, (Subscript[n, p]) of which form a product P, as well as (Subscript[n, e]) enzymes, (Subscript[n, c]) competitive inhibitors and (Subscript[n, u]) uncompetitive inhibitors, the rate of product formation is

\(TraditionalForm\`d[P]\/dt = ∑\+\(z = 1\)\%\(n\_p\)\(∑\+\(a = 1\)\%\(n\_e\)\(\(V\_\(max\_\(E\_a\)\)\) \([S\_z]\)\/K\_\(S\_z, \ E\_a\)\)\/\(1 + ∑\+\(b = 1\)\%\(n\_s\)\([S\_b]\)\/K\_\(S\_b, \ E\_a\) + ∑\+\(c = 1\)\%\(n\_c\)\([I\_c]\)\/K\_\(I\_c, \ E\_a\) + ∑\+\(u = 1\)\%\(n\_u\)\(∑\+\(b = 1\)\%\(n\_s\)\(\([S\_b]\)\/K\_\(S\_b, \ E\_a\)\) \([I\_u]\)\/K\_\(I\_u, \ \(E\_a\) S\_b\)\)\)\)\)

We are still trying to incorporate the noncompetitive rate law into the the general form.

Nontrivial Michaelis-Menten example

The neurotoxic amyloid protein Aβ that accumulates in the brains of Alzheimer's disease patients is produced from its precursor protein APP by the following mechanism.

Several drug companies are pursuing inhibitors of γ-secretase as a therapy for Alzheimer's disease. Here is how we could model the biochemistry for this part of the disease process.

\[KernelIcon]

model={
 (Overscript[↦ , (Subscript[V, APP])])APP,
APP(Overscript[→, MM[αSec,(Subscript[K, m, α]),(Subscript[V, max, α])]])sAPPα + 
                C83,
APP(Overscript[→, MM[βSec,(Subscript[K, m, β]),(Subscript[V, max, β])]])sAPPβ + 
                C99,
C83 (Overscript[→ , MM[γSec, 
                    (Subscript[K, m, γ, p3]), 
                    (Subscript[V, max, γ, p3])]])p3 + AICD,
C99 (Overscript[→ , MM[γSec, 
                    (Subscript[K, m, γ, Aβ40]), 
                    (Subscript[V, max, γ, Aβ40])]])Aβ40+AICD,
C99 (Overscript[→ , MM[γSec, 
                    (Subscript[K, m, γ, Aβ42]), 
                    (Subscript[V, max, γ, Aβ42])]])Aβ42+AICD,
inhib (Overscript[⇥ , MMcomp[γSec, (Subscript[K, m, inhib])]])};
{vars, {odes, ics}} = 
        KineticModel[model, t];
Cases[
          odes, (([Aβ40 | Aβ42])^′)[t] == _] // 
        ToColumns[#, 1] & // 
      Grid[#, RowSpacings -> 3, 
          ColumnAlignments -> "=="] &

Pharmacokinetics

\[KernelIcon]

SetOptions[StandardizeUnits, StandardAmount -> Nano Mole, 
      StandardTime -> Hour]

\[KernelIcon]

model = {
⟨ (Overscript[↦ , drugAdmin])drug⟩_gut, 
⟨drug⟩_gut (Overscript[⇒ , (Subscript[k, in])])⟨drug⟩_plasma, 
⟨drug (Overscript[→ , (Subscript[k, m])])metabolite⟩_plasma, 
⟨drug (Overscript[⇥ , (Subscript[k, elim, drug])])⟩_plasma, 
⟨metabolite (Overscript[⇥ , (Subscript[k, elim, met])])⟩_plasma
};
{vars, {odes, ics}} = KineticModel[model, t];
odes // ToColumns[#, 1] & // 
      Grid[#, RowSpacings -> 3, 
          ColumnAlignments -> "=="] &

\[KernelIcon]

rateRules = {(Subscript[k, in]) -> 0.08  (Hour^(-1)), 
            (Subscript[k, m]) -> 0.15  (Hour^(-1)), 
            (Subscript[k, elim, drug]) -> 0.04  (Hour^(-1)), 
            (Subscript[k, elim, met]) -> 0.02475  (Hour^(-1))} // 
          StandardizeUnits;

Single dose

\[KernelIcon]

?Bolus

\[KernelIcon]

adminFunc = Bolus[0.005, 8,  t]

\[KernelIcon]

<< DifferentialEquations`Impulse`

\[KernelIcon]

solnSingle = 
        NDSolve[{odes, ics} /. drugAdmin -> adminFunc /. rateRules // 
              ShedUnits, vars, {t, 0, 600}, Method -> Impulse] // 
          First;
plotPK[adminFunc, 
        solnSingle, "Single Dose on Day 1"];

Multiple doses

\[KernelIcon]

adminFunc = (∑_(i = 1)^28) Bolus[0.005, 8 + (i - 1) 12,  
            t];
solnRepeated = 
        NDSolve[{odes, ics} /. drugAdmin -> adminFunc /. rateRules // 
              ShedUnits, vars, {t, 0, 28×24}, Method -> Impulse] // 
          First;
plotPK[adminFunc, 
        solnRepeated, {28, 0.015}, "BID for 14 days"];

Multiple dose, non-compliant patient

Completely forgot morning dose on day 4; missed evening dose on day 8 and took 2 doses morning of day 9.

\[KernelIcon]

adminFunc = ((∑_(i = 1)^28) Bolus[0.005, 
                8 + (i - 1) 12,  t]) - 
          Bolus[0.005, 8 + (8 - 1) 12, t] - 
          Bolus[0.005, 8 + (17 - 1) 12, t] + 
          Bolus[0.005, 8 + (18 - 1) 12, t];
solnMissed = 
        NDSolve[{odes, ics} /. drugAdmin -> adminFunc /. rateRules // 
              ShedUnits, vars, {t, 0, 28×24}, Method -> Impulse] // 
          First;
plotPK[adminFunc, 
        solnMissed, {28, 
          0.015}, "BID for 14 days with missed and doubled \
doses"];

Bolus + Infusion

\[KernelIcon]

?Infusion

\[KernelIcon]

adminFunc = Bolus[0.006, 8, t] + Infusion[0.02, 8, 56, t]

\[KernelIcon]

solnLoadingDose = 
        NDSolve[{odes, ics} /. drugAdmin -> adminFunc /. rateRules // 
              ShedUnits, vars, {t, 0, 12×24}, MaxStepSize -> 1, 
            Method -> Impulse] // First;
plotPK[adminFunc, 
        solnLoadingDose, {12, 0.015}, "Bolus + Infusion"];

Parameter optimization

Very often we want to optimize the rate constants against a set of experimental data. For example, flumazenil (FMZ), a specific and competitive antagonist at the central benzodiazepine receptor, is used as a radiotracer in PET imaging studies of the brain. It is metabolized in the body to a carboxylic acid (COOH), and understanding the rate at which this happens is important.

The data we will fit against came from Debruyne, D., P. Abadie, et al. (1991). “Plasma pharmacokinetics and metabolism of the benzodiazepine antagonist (11C)Ro 15-1788 (flumazenil) in baboon and human during positron emission tomography studies.” Eur. J. Drug Metab. Pharmacokinet. 16(2): 141-52.

\[KernelIcon]

With[{data = << baboonData.m}, colNames = First[data]; 
        baboonData = Rest[data]];
Grid[Prepend[baboonData, colNames], 
      ColumnAlignments -> Decimal]
plotFit1["Baboon"];

\[KernelIcon]

\(\(model = {\(\( ↦ \&\(Bolus[100, 0, 
                  t]\)\)\(⟨FMZ⟩\_1\)\), 
⟨FMZ⟩\_1\( ⇔ \+\(K\_\(2, 1\)\)\%\(K\_\(1, 
                    2\)\)\)⟨FMZ⟩\_2, 
\(\(⟨FMZ⟩\_1\)\( ⇥ \&\(K\_e\)\)\), 
⟨FMZ\( → \&\(K\_m\)\)COOH⟩\_1, 
\(\(⟨COOH⟩\_1\)\( ⇥ \&\(K\_\(m, e\)\)\)\)};\)\)
{vars, {odes, ics}} = 
        KineticModel[model, t, IgnoreVolumes -> False, 
          InitialTime -> -5];
odes // List // Thread // 
      Grid[#, RowSpacings -> 4, ColumnAlignments -> Left, 
          ColumnWidths -> 50] &

Here is an initial guess at the parameters:

\[KernelIcon]

baboonParams = {injectedDose -> 
            12.95 Milli Curie×(((10^(-6)) Mole)/(388.5 Milli Curie)), (Subscript[K, e]) -> 0.086 (Minute^(-1)), 
          (Subscript[K, 1, 2]) -> 0.205 (Minute^(-1)), 
          (Subscript[K, 2, 1]) -> 0.112 (Minute^(-1)), 
          (Subscript[K, m]) -> 0.022 (Minute^(-1)), 
          (Subscript[K, m, e]) -> 0.375 (Minute^(-1)), 
          CompartmentVolume[1] -> 2.2 Liter, 
          CompartmentVolume[2] -> 3.8 Liter};
soln = 
        NDSolve[{odes, ics} /. baboonParams // ShedUnits, 
            vars,  {t,  -5,  60}, Method -> Impulse] // 
          First;
plotFit2[soln, "Baboon"];

Data fitting

There are a number of things that need to be addressed before we can do the fit.

({{1., Problem:, \(As\ biologists\ we\ often\ want\ to\ compare\ the\ fitted\ parameters\ from\ one\ experiment\ to\ those\ from\ another\ and\ make\ a\ statistical\ test\ of\ significance\ on\ their\ \(\(difference\)\(.\)\)\)}, { , Solution:, To get confidence intervals of the fitted parameters we will use NonlinearRegress.}, {2., Problem:, We have two response variables:[FMZ][t] and[COOH][t].}, { , Solution:, We will recast the data table to have 1 response variable column and add an indicator variable column to specify to which response an observation corresponds .  The response expression will be a Piecewise function based on the indicator variable.}, {3., Problem:, We don't have an analytical expression for the response variables in terms of the adjuatable parameters (NDSolve returns InterpolatingFunctions).}, { , Solution:, We will take derivatives of the ODEs with respect to the parameters and teach Mathematica how to use them  with the derivative function D.}})

Derivatives with respect to parameters

We have generalized the method of Carl Woll (in Paul Abbott, "In and Out" The Mathematica Journal, 2004, 9, 310-313), and implemented it in the package ODEDerivatives.m

\[KernelIcon]

<< Chemistry`ODEDerivatives`

\[KernelIcon]

?DefineODEDerivatives

Adjustable and fixed parameters

\[KernelIcon]

adjustableParams = 
      Cases[baboonParams, 
        HoldPattern[
          (Subscript[K, e]) | (Subscript[K, m]) | (Subscript[K, m, e]) | (Subscript[K, 1, 2]) | 
              (Subscript[K, 2, 1]) -> _]]
adjustableParamNames = First /@ adjustableParams;

\[KernelIcon]

fixedParams = 
      Cases[baboonParams, 
        HoldPattern[CompartmentVolume[_] -> _ // Evaluate]]
fixedParamNames = First /@ fixedParams;

\[KernelIcon]

DefineODEDerivatives[vars,{odes, ics},{t, -5, tmax},ℱ,adjustableParamNames,fixedParamNames,Method -> Impulse]

Note that the last argument allows one to supply options to NDSolve.

The rules for the solution of the ODEs and its derivatives with respect to the parameters are stored in the SubValues of Derivative and ℱ. Note the use of dynamic programming to make the computation more efficient.

\[KernelIcon]

SubValues[
        Derivative][[3]]

\[KernelIcon]

SubValues[ℱ][[1]]

Set up the fitting data table and model

The third and fourth columns of the data table have the response variables. We set the indicator variable i to 1 for [FMZ][t] and to 2 for [COOH][t].

\[KernelIcon]

fitData=Join@@({{#1,1,#3},{#1,2,#4}}&)@@@baboonData

Here is the response function using Piecewise

\[KernelIcon]

fitFunc=((Piecewise[{{ℱ[[⟨FMZ⟩_1]][(Subscript[K, e]), (Subscript[K, 1, 2]), (Subscript[K, 2, 1]), (Subscript[K, m]), 
                        (Subscript[K, m, e]), 
CompartmentVolume[1], 
                        CompartmentVolume[2], 60][t], i == 1}, {ℱ[[⟨COOH⟩_1]][(Subscript[K, e]), (Subscript[K, 1, 2]), (Subscript[K, 2, 1]), (Subscript[K, m]), 
                        (Subscript[K, m, e]), 
CompartmentVolume[1], 
                        CompartmentVolume[2], 60][t], i == 2}}]))/.fixedParams//ShedUnits

Fitting the data

\[KernelIcon]

<< Statistics`NonlinearFit`

We have added two new RegressionReport options to Regress and NonlinearRegress: SinglePredictionCI and MeanPredictionCI.

\[KernelIcon]

?SinglePredictionCI

\[KernelIcon]

fit = 
        NonlinearRegress[fitData, fitFunc, {t, i}, 
          List @@@ adjustableParams // ShedUnits, 
          RegressionReport -> {SummaryReport, BestFit, 
              SinglePredictionCI}];
ANOVATable /. fit
ParameterCITable /. fit
plotFit3[soln, "Baboon"];

Here is how we can display the confidence interval about the fitted curves:

\[KernelIcon]

f[i_] = BestFit /. fit
{lci[i_], uci[i_]} = SinglePredictionCI /. fit;

\[KernelIcon]

DisplayTogether[ListPlot[Cases[fitData, {t_, 1, y_} :> {t, y}], 
            PlotStyle -> {Red, AbsolutePointSize[5]}, 
            Axes -> False, Frame -> True, 
            FrameLabel -> {StyleForm["Time (min)", 
                  FontFamily -> "Times", FontSize -> 12],  
                StyleForm["% injected dose", FontFamily -> "Times", 
                  FontSize -> 12]}, 
            PlotRange -> {{-2, 62}, {-2, 45}}, 
            PlotLabel -> 
              StyleForm["Baboon", FontFamily -> "Times", 
                FontSize -> 16]],ListPlot[Cases[fitData, {t_, 2, y_} :> {t, y}], 
            PlotStyle -> {Blue, AbsolutePointSize[5]}],Plot[{f[1], lci[1], uci[1], f[2], lci[2], 
                uci[2]}, {t, -5, 60}, 
              PlotStyle -> {{Red}, {Red, Dashed}, {Red, 
                    Dashed}, {Blue}, {Blue, Dashed}, {Blue, Dashed}}],\(Graphics[{MapThread[{Text[
                      StyleForm[#1, FontFamily -> "Times", 
                        FontSize -> 12], {45, 
                        30 - 4 \((#3 - 1)\)}, {\(-1\), 0}], Sequence @@ #2, 
                    Line[{{36, 30 - 4 \((#3 - 1)\)}, {43, 
                          30 - 4 \((#3 - 
                                  1)\)}}]} &, {{\([⟨FMZ⟩\_1]\), \*""[⟨FMZ⟩_1] 95% CI"", \([⟨COOH⟩\_1]\), \*""[⟨COOH⟩_1] 95% CI""}, {{Red}, {Red, 
                      Dashed}, {Blue}, {Blue, Dashed}}, Range[4]}]}]\),ImageSize -> 7×72];

Other applications

The tools in ChemicalKinetics.m can be used wherever one is working with compartmental models. Classic examples include:

Others' work

We have not been the first, and I hope we are not the last, to use Mathematica to develop kinetic modeling tools. Here is a list of publications and products:

We apologize for any omissions, and would appreciate being made aware of them.

Acknowledgements

Several people have contributed ideas and time over the past few years to the work described here. I would like to especially thank

Initialization

Off[General::spell1, General::spell]
ToFileName[{$HomeDirectory, "Presentations", "2006", "WRI \
Technology Conference"}] // SetDirectory;
<< Graphics`
<< Chemistry`ChemicalKinetics`
<< Statistics`NonlinearFit`
DarkGray = GrayLevel[1/3] // N;

Formatting

ToColumns[list_List, nCols_Integer?Positive, fill_:  Null] := 
      With[{n = Ceiling[Length[list]/nCols]}, 
        Partition[list, n, n, {1, 1}, {fill}] // 
          Transpose]
ToRows[list_List, nRows_Integer?Positive, fill_:  Null] := 
      With[{n = Ceiling[Length[list]/nRows]}, 
        Partition[list, n, n, {1, 1}, {fill}]]

ClearAll[MakeBoxes]

\(MakeBoxes[
        l : {H2 + N2O2\( → \&\(k\_1\)\)H2O + N2O, 
            2\ NO\( ⇌ \+\(k\_5\)\%\(k\_4\)\)N2O2, 
            H2 + N2O\( → \&\(k\_6\)\)H2O + N2}, StandardForm] := 
      StyleBox[listBox[l], FontColor -> Red]\)

\(MakeBoxes[
        l : {\(HCO3\^-\) + \(Na\^+\)\( ⇌ \+\(k\_3\)\%\(k\_2\)\)NaHCO3, 
            NaOH\( ⇌ \+\(k\_8\)\%\(k\_7\)\)\(OH\^-\) + \(Na\^+\), 
            CO2 + \(OH\^-\)\( ⇌ \+\(k\_10\)\%\(k\_9\)\)\(HCO3\^-\)}, StandardForm] := StyleBox[listBox[l], FontColor -> Blue]\)

insertCommas[l_List] := 
      Fold[Join[#1, {",", #2}] &, {First[l]}, Rest[l]]

listBox[l_List] := 
      RowBox[{"{", 
          RowBox[MakeBoxes[#, StandardForm] & /@ l // 
              insertCommas], "}"}]

MakeBoxes[x:Times[Power[Liter, -2], Power[Mole, 2], 
                Power[Second, -1]],StandardForm]:=StyleBox[fractionBox[x], FontColor -> Red]

fractionBox[x_] := 
      With[{n = Numerator[x], d = Denominator[x]}, 
        FractionBox[MakeBoxes[n, StandardForm], 
          MakeBoxes[d, StandardForm]]]

MakeBoxes[r:(\(k\_2 -> 5.`\/Second\)),StandardForm]:=StyleBox[ruleBox[r], FontColor -> Blue]

ruleBox[lhs_ -> rhs_] := 
      RowBox[{MakeBoxes[lhs, StandardForm], "->", 
          MakeBoxes[rhs, StandardForm]}]

Plotting

Clear[plotModel]
\(\(plotModel[soln_] := 
        Plot[
{
\([NO]\)[
                  t], 
\([H2]\)[
                  t], 
\([N2O2]\)[
                  t], 
\([N2O]\)[
                  t], 
\([N2]\)[
                  t], 
\([H2O]\)[t]
} /. 
              soln // Evaluate, 
{t, 0, 
            600}, 

ImageSize -> 
            600, 
Axes -> 
            None, 
Frame -> True, 
          FrameLabel -> {StyleForm["Time (min)", 
                FontFamily -> "Times", FontSize -> 12], \ 
              StyleForm["Concentration (M)", FontFamily -> "Times", 
                FontSize -> 
                  12]}, 
FrameTicks -> {Automatic, 
              Automatic}, 

PlotLabel -> StyleForm[\*""2NO + 2(Subscript[H, 2]) → (Subscript[N, 2]) + 2(Subscript[H, 2])O"", FontFamily -> "Times", 
              FontSize -> 
                16], 
PlotRange -> {Automatic, {\(-0.0002\), 0.025}}, 
PlotPoints -> 
            500, 
PlotStyle -> {{RGBColor[1\/2, 0, 
                  1\/2], Solid}, {Black, Solid}, {RGBColor[1, 0, 1], 
                Dashed}, {RGBColor[1\/3, 0, 2\/3], 
                Dashed}, {RGBColor[0, 0, 1], Solid}, {Gray, 
                Solid}}, 

PlotLegend -> {
StyleForm["[NO]", FontFamily -> "Times", 
                FontSize -> 12], 
StyleForm[\([H2]\), 
                FontFamily -> "Times", 
                FontSize -> 
                  12], 
StyleForm[\([N2O2]\), 
                FontFamily -> "Times", 
                FontSize -> 12], 
StyleForm[\([N2O]\), 
                FontFamily -> "Times", 
                FontSize -> 12], 
StyleForm[\([N2]\), 
                FontFamily -> "Times", 
                FontSize -> 12], 
StyleForm[\([H2O]\), 
                FontFamily -> "Times", FontSize -> 12]}, 
          
          LegendPosition -> {0.25, \(-0.2\)}, 
LegendShadow -> 
            None, 
LegendSize -> {1, 
              0.5}, 
LegendSpacing -> 
            0.001
];\)\)

Clear[plotPK]
\(\(plotPK[drugAdmin_, soln_, xymax_: {12, 0.006}, title_String] := 
        With[{vars = {\([⟨drug⟩\_gut]\), \([⟨drug⟩\_plasma]\), \([⟨metabolite⟩\_plasma]\)}, 
            xmax = 24  First[xymax], ymax = Last[xymax], 
            colors = {ForestGreen, Red, Blue}}, 
          DisplayTogether[
            Plot[
\((vars[t] // Through)\) /. soln // 
                Evaluate, 
{t, 0, 
                xmax}, 

Axes -> 
                None, 
Frame -> True, 
              FrameLabel -> {StyleForm["Days Elapsed", 
                    FontFamily -> "Times", FontSize -> 12], \ 
                  StyleForm["Concentration", FontFamily -> "Times", 
                    FontSize -> 
                      12]}, 
FrameTicks -> {Table[{i, 
                      If[Mod[i, 2\ 24] == 0, i\/24, ""]}, {i, 0, 
                      xmax, 24\ }], Automatic, 
                  Table[{i, ""}, {i, 0, xmax, 24\ }], 
                  Automatic}, 

PlotLabel -> StyleForm[title, FontFamily -> "Times", 
                  FontSize -> 16], 
PlotPoints -> 
                50, 
PlotStyle -> 
                colors
], 

With[{dmax = 
                  Join[Cases[
                        drugAdmin // 
                          Expand, \((d_\ DiracDelta[_])\) :> 
                          d/\((xmax\/500)\), {0, ∞}], 
                      Cases[drugAdmin // 
                          Expand, \((d_\ UnitStep[_])\) :> 
                          d, {0, ∞}]] // Max}, 
              Plot[\((0.6 \( ymax\/\(4  dmax\)\) drugAdmin - 0.8  ymax\/4)\) /. 
                    d_\ DiracDelta[t + x_] :> 
                      d/\((xmax\/500)\) \((UnitStep[t + x] - 
                            UnitStep[t + x - xmax\/500])\) // Evaluate, {t, 
                  0, xmax}, PlotPoints -> 750, 
                PlotStyle -> {DarkGray}]], 

Graphics[{MapThread[{Text[
                        StyleForm[#1, FontFamily -> "Times", 
                          FontSize -> 12], {0.75  xmax, 
                          0.8  ymax - 0.075  ymax \((#3 - 1)\)}, {\(-1\), 
                          0}], #2, 
                      Line[{{0.63  xmax, 
                            0.8  ymax - 0.075  ymax \((#3 - 1)\)}, {0.73  
                              xmax, 0.8  ymax - 
                              0.075  ymax \((#3 - 1)\)}}]} &, {vars, colors, 
                    Range[Length[vars]]}], 
                Text[StyleForm["drugAdmin", FontFamily -> "Times", 
                    FontSize -> 12], {0.75  xmax, \(-0.6\) 
                      ymax\/4}, {\(-1\), 
                    0}]}], 

ImageSize -> 600, 
PlotRange -> {Automatic, {\(-\(ymax\/4\)\), ymax}}
]];\)\)

Clear[plotFit1]
\(\(plotFit1[xymax_: {62, 45}, title_String] := 
        With[{vars = {\([⟨FMZ⟩\_1]\), \([⟨COOH⟩\_1]\), \([⟨FMZ⟩\_1]\) + \([⟨COOH⟩\_1]\)}, 
            colors = {Red, Blue, Black}}, 
          DisplayTogether[
ListPlot[
              baboonData[[All, {1, 3}]], 
              PlotStyle -> {colors[[1]], AbsolutePointSize[5]}, Axes -> False, 
              Frame -> True, 
              FrameLabel -> {StyleForm["Time (min)", 
                    FontFamily -> "Times", FontSize -> 12], \ 
                  StyleForm["% injected dose", FontFamily -> "Times", 
                    FontSize -> 12]}, 
              PlotRange -> {{\(-2\), 62}, {\(-2\), 45}}, 
              PlotLabel -> 
                StyleForm[title, FontFamily -> "Times", 
                  FontSize -> 16]], 
ListPlot[
              baboonData[[All, {1, 3}]], 
              PlotStyle -> {colors[[1]], AbsolutePointSize[5]}, 
              PlotJoined -> 
                True], 

ListPlot[
              baboonData[[All, {1, 4}]], 
              PlotStyle -> {colors[[2]], AbsolutePointSize[5]}], 
ListPlot[
              baboonData[[All, {1, 4}]], 
              PlotStyle -> {colors[[2]], AbsolutePointSize[5]}, 
              PlotJoined -> 
                True], 

ListPlot[
              baboonData[[All, {1, 2}]], 
              PlotStyle -> {colors[[3]], AbsolutePointSize[5]}], 
ListPlot[
              baboonData[[All, {1, 2}]], 
              PlotStyle -> {colors[[3]], AbsolutePointSize[5]}, 
              PlotJoined -> 
                True], 

Graphics[{MapThread[{Text[StyleForm[#1, FontFamily -> "Times", 
                          FontSize -> 12], {43, 
                          30 - 4 \((#3 - 1)\)}, {\(-1\), 0}], #2, 
                      Line[{{34, 30 - 4 \((#3 - 1)\)}, {41, 
                            30 - 4 \((#3 - 1)\)}}]} &, {vars, colors, 
                    Range[Length[
                        vars]]}]}], 

ImageSize -> 
              7×72]];\)\)

Clear[plotFit2]
\(\(plotFit2[soln_, xymax_: {62, 45}, title_String] := 
        With[{vars = {\([⟨FMZ⟩\_1]\), \([⟨COOH⟩\_1]\), \([⟨FMZ⟩\_2]\)}, 
            colors = {Red, Blue, DarkGreen}}, 
          DisplayTogether[
ListPlot[
              baboonData[[All, {1, 3}]], 
              PlotStyle -> {colors[[1]], AbsolutePointSize[5]}, Axes -> False, 
              Frame -> True, 
              FrameLabel -> {StyleForm["Time (min)", 
                    FontFamily -> "Times", FontSize -> 12], \ 
                  StyleForm["% injected dose", FontFamily -> "Times", 
                    FontSize -> 12]}, 
              PlotRange -> {{\(-2\), 62}, {\(-2\), 45}}, 
              PlotLabel -> 
                StyleForm[title, FontFamily -> "Times", 
                  FontSize -> 
                    16]], 

ListPlot[
              baboonData[[All, {1, 4}]], 
              PlotStyle -> {colors[[2]], 
                  AbsolutePointSize[
                    5]}], 

Switch[
              soln, {\((_ -> _InterpolatingFunction)\) .. }, 
              Plot[Through[vars[t]] /. soln // Evaluate, {t, \(-5\), 60}, 
                PlotStyle -> colors], _, 
              Graphics[{}]], 

Graphics[{MapThread[{Text[
                        StyleForm[#1, FontFamily -> "Times", 
                          FontSize -> 12], {47, 
                          30 - 4 \((#3 - 1)\)}, {\(-1\), 0}], #2, 
                      Line[{{38, 30 - 4 \((#3 - 1)\)}, {45, 
                            30 - 4 \((#3 - 1)\)}}]} &, {vars, colors, 
                    Range[Length[
                        vars]]}]}], 

ImageSize -> 
              7×72]];\)\)

Clear[plotFit3]
\(\(plotFit3[soln_, xymax_: {62, 45}, title_String] := 
        With[{vars = {\([⟨FMZ⟩\_1]\), \([⟨COOH⟩\_1]\), \([⟨FMZ⟩\_2]\)}, 
            colors = {Red, Blue, DarkGreen}}, 
          DisplayTogether[
ListPlot[
              Cases[fitData, {t_, 1, y_} :> {t, y}], 
              PlotStyle -> {colors[[1]], AbsolutePointSize[5]}, Axes -> False, 
              Frame -> True, 
              FrameLabel -> {StyleForm["Time (min)", 
                    FontFamily -> "Times", FontSize -> 12], \ 
                  StyleForm["% injected dose", FontFamily -> "Times", 
                    FontSize -> 12]}, 
              PlotRange -> {{\(-2\), 62}, {\(-2\), 45}}, 
              PlotLabel -> 
                StyleForm[title, FontFamily -> "Times", 
                  FontSize -> 
                    16]], 

ListPlot[
              Cases[fitData, {t_, 2, y_} :> {t, y}], 
              PlotStyle -> {colors[[2]], 
                  AbsolutePointSize[
                    5]}], 

Plot[\(\((\(\(\(ℱ[#]\)[K\_e, K\_\(1, 2\), K\_\(2, 1\), K\_m, K\_\(m, e\), 
                              CompartmentVolume[1], CompartmentVolume[2], 
                              60]\)[t] &\) /@ 
                        vars)\) /. \((BestFitParameters /. fit)\)\) /. \((fixedParams // ShedUnits)\) // 
                Evaluate, {t, \(-5\), 60}, 
              PlotStyle -> 
                colors], 

Graphics[{MapThread[{Text[
                        StyleForm[#1, FontFamily -> "Times", 
                          FontSize -> 12], {47, 
                          30 - 4 \((#3 - 1)\)}, {\(-1\), 0}], #2, 
                      Line[{{38, 30 - 4 \((#3 - 1)\)}, {45, 
                            30 - 4 \((#3 - 1)\)}}]} &, {vars, colors, 
                    Range[Length[
                        vars]]}]}], 

ImageSize -> 
              7×72]];\)\)

Reset $Line

$Line = 0